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FOREWORD 


This report is submitted by Rocketdyne to the National Aeronautics and Space 
Administration Office of Advanced Research and Technology as the final 
technical report of the research performed under Contract NAS 7 -470 during 
the period 13 June 1966 through 12 May 1968 . Technical manager of this 
contract is Dr. R. H. Wilson, Chief, Applied Mathematics Branch, NASA 
Headquarters, Washington, D.C. 


ABSTRACT 

A multistep predictor-corrector method for the numerical solution of ordinary 
differential equations is developed. The difference equations employed are 
generalizations, for the case of variable mesh spacing, of previous formulas 
requiring fixed step size. In addition to retaining the high local accuracy 
and convergence properties of the earlier methods, the variable mesh method 
is developed in a form conducive to the generation of effective criteria for 
the selection of subsequent step sizes in the step by step solution of differ- 
ential equations . These criteria are based on considerations of truncation 
error, convergence of corrector iterations, and an extensive treatment of 
relative numerical stability. The algorithm has been tested extensively and 
compared with other methods. The results of the comparison favor the new 
method. 


This report also discusses an extension of the variable mesh multistep method 
for the case of stiff equations and application of the variable mesh approach 
to partial differential equations. 
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VARIABLE MESH METHODS FOR 


DIFFERENTIAL EQUATIONS 

1. INTRODUCTION. A great deal of research effort has been directed 
toward the numerical solution of first order nonlinear ordinary differen- 
tial equations because of the practical importance of such problems . The 
most widely used numerical methods that have been developed for these pro- 
blems provide approximate values of the solution at discrete points accord- 
ing to a stepwise computation beginning at an initial point for which the 
solution is known. These methods are called one-step methods if the cal- 
culation of the solution at a given point depends explicitly on values of 
the solution and one or more of its derivatives at only one previous point. 
Multistep methods require values at two or more previous points. One-step 
(Runge-Kutta) methods are very convenient because the step increments can 
be changed readily from step to step as desired and because the solution 
in the initial steps is calculated with the same formulas as used in sub- 
sequent steps. Multistep methods, although less convenient, are usually 
more efficient because, by making use of the calculations of more than one 
previous step, less computer time is required to achieve the same accuracy 
as achieved with a one -step method. 

The research reported here was directed toward the development and testing 
of variable mesh multistep methods which not only preserve the efficiency 
due to the multistep structure but improve this efficiency by permitting as 
much freedom in the variation of the step increments as is afforded by one- 
step methods . Care was taken to formulate the basic difference equations 
in a manner conducive to the development of effective criteria for selecting 
the variable mesh increments as the calculation progresses . In the follow- 
ing pages the basic algorithm is described and the analysis and practical 
considerations justifying the mesh criteria are presented. The mesh criteria 
were subjected to extensive numerical testing, and in addition, the algorithm 
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was compared with known methods in the numerical solution of selected dif- 
ferential equations . The results of this and other experimental work are 
Busmarized in later sections. 

The problem of starting the computation, that is, the requirement of com- 
puting the solution at the first few points by a separate technique in 
order to initialize multistep methods, is not emphasized here for two 
reasons . First, because of the variable mesh formulation, the calculation 
is only initialized once and never has to be restarted as would be required 
in changing the step size while using a fixed step size, multistep method. 
In the second place, fairly general starting procedures are readily avail- 
able for incorporation with the variable mesh method because the step in- 
crements used in the starting procedure can be smaller than those used in 
the subsequent calculations . For example, the starting procedure outlined 
in Cl] for the variable mesh method consists of simply using the one-step 
Adams -Bashf orth/Adams -Moulton formulas for the first step, the two-step 
formulas for the second step, and the three-step formulas for third step. 
Die same step size is used for each of these three initial steps, and it 
is chosen small enough to yield the desired accuracy at the first point. 
There is little danger of exceeding this error at the second and third 
points since higher order formulas are used. 

In addition to the basic variable mesh multistep method for ordinary dif- 
ferential equations, this report discusses an extension of the method for 
stiff equations and application of the variable mesh approach to partial 
differential equations. 
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2. VARIABLE MESH MULTISTEP FORMULAS. The Initial value problems 
considered are represented by differential equations of the form 

g-F(x,y) , (1) 

with initial condition y(x Q ) = y Q - Equation (l) represents a single 
differential equation; however some special considerations required 
for systems of differential equations will also be given in the sequel. 

It is assumed at the outset that F is continuous and satisfies the 
Lipschitz condition that guarantees the existence of a unique, continu- 
ous and differentiable solution [2, p. 15]. The continuity of higher 
derivatives will be required later in the discussion of truncation 
error . 

We will use the usual notation in which y denotes the computed value 

of y(x ) and y f denotes F(x , y ) • It is assumed that the computed 
n n n' n 

solution is obtained recursively by one or more formulas of the following 
type: 


r»l " a o y n + Vn-l + a 2 y n-2 + a 3 y n-3 


+ h ( b -l y n + l + Vn + b l y A-l + Vi-2 + b 3 y n-3 } * 


( 2 ) 


Here h denotes the current step size, x . - x , and is permitted to 

n+1 n 

vary with n. The coefficients a^ and b^ are also variable and it 
will be convenient later to express them in terms of mesh parameters ct, 
B, and v defined by 


° - ( Wi )/h ' 
e " (x n“ x n-2 )/h ' 

Y * < x n“ x n-3 )/h * 


( 5 ) 
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We require y > P > « 


For the case of fixed step size, a, 0, and y have the constant values 
1, 2, and 5, respectively. Henrici [2, pp. 218, 224] has defined condi- 
tions of consistency and stability for fixed h and has shown that they 
are necessary and also sufficient, when taken together, for convergence 
of y & to y(x Q ) as h-*0. The stability condition requires that no 
root of the equation 

P* ~ a o P 5 - ajP 2 - a 2 p - a ? = 0 (4) 

exceed one in modulus and roots of unit modulus must be simple . The 
consistency condition requires that equation (2) be exact if y(x) is 
either constant or linear. 

An analogous equivalence theorem holds in the variable mesh case. It 
follows Immediately that the stability and consistency conditions are 
necessary for convergence in the variable mesh case because they are 
necessary in the special case of 'fixed mesh. Henrici* s proof of suffi- 
ciency has been generalized by the author to account for the case of 
variable step size, but is emitted here because of its length. For this 
case the two conditions are required to hold for each different step size 
used in the integration. 

The validity of the generalized equivalence theorem is not restricted to 
difference equations with only the number of terms actually shown in (2). 
However since we will restrict the present discussion to fourth order 
methods— that is, methods with error terms proportional to the fifth power 
of h— the terms 6hown are adequate. The optimum order to use in a given 
application depends heavily on the degree of accuracy desired, but fourth 
order is a reasonable c 005) remise for medium range accuracy— say two to 
six significant figures. With fixed step size it is often desirable to 
vary the order within a given application in order to maintain a desired 
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accuracy. The variable mesh procedure, on the other hand, has the ad- 
vantage of achieving the same objective without switching from formulas 
of one order to those of another. 


Both explicit (predictor) and implicit (corrector) variable mesh formulas 
are used. The explicit equation has b_^ = 0 while the implicit usually 
can be solved by iteration. The coefficients in ( 2 ) for the two formulas 
are determined in part by requiring satisfaction of the stability and 
consistency conditions. By requiring exactness for F(x,y) = 0, the rela- 
tion, a Q + a^ + a 2 + a^ = 1, is imposed, from which it follows that one 
root of equation (k) is unity. The other three roots ("parasitic " \ arising 
because a fourth order difference equation is used in place of a first 
order differential equation, are all zero if we select a Q = 1 , a^ = ~ 



The consistency condition is satisfied by the additional requirement 
of exactness for y = x ” x n > which yields, for the predictor, b Q + b^ + 
b^ + b^ = 1 . Combining this with the requirements of exactness for 

$( x ~ x a ) 2 > 3^ x ~ x n ^’ £( x-x n )^> the b i of the P redictor > 


y n+l = y n + h <Vn + Vn-1 + Vn-2 + h 3 y n-3 } ’ 


(5) 


are determined recursively as follows: 

2(2+3a) (fi+oc) + 3(l-2a 2 ) 
b 3 = i2y(y-a)(/3-y) 


2+3 a-6j/(y-a)bj 

b 2 = 6 

b l = - ^( 1+ 2yb 3 +2^h 2 ) 


b 0 - 1 - b 3 - b 2 - b l 


( 6 ) 
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Similarly, a corrector of the form 


y n+l = y n + h < d -l y n+l + d o y n + d l y n-l + d 2 y n-2> 
is found with coefficients 



For the special case of fixed step size, the above predictor and corrector 
formulas reduce to the widely accepted Adams -Bashforth and Adams-Moulton 
formulas respectively. In this connection one is reminded of the formulas 
presented by Nordsieck in a paper which, like the present paper, also em- 
phasizes the advantages of changing step size [3]- Although the algorithm 
of Nordsieck is substantially different from that presented here, it is 
similar in the sense that his basic integration formulas are equivalent 
to the Adams formulas. However, the formulation used by Nordsieck appears 
to be much less conducive to the development of effective mesh selection 
criteria than is the formulation presented above. This claim is corroborated 
by evidence obtained when both methods, complete with their respective 
recommended mesh selection criteria, were applied to selected differential 
equations. (This work is described in more detail in a later section.) 

This deficiency of the Nordsieck method may result from his assumption 
of "fixed point" operations rather than the more commonly used "floating 
point, " as also suggested by Levis and Stovall in a paper which appeared 
too late for incorporation in the work reported here [4]* 


6 



41 


Assuming continuous higher derivatives of F(x, y), it is evident upon com- 
paring equation ( 5 ) with an appropriate Taylor Series representation for 
y*(x that the truncation error in ( 5 ) can be represented as 



0(h 6 ), 


where the coefficient depends on a, and y. If we consider the 
residual error resulting from the application of ( 5 ) to the polynomial 
(x-x )^, we find that 

p n = 1 - 5 (b ia 4 +b 2 /3 4 +b 3 y 4 ) 

= 1+12 [3 (a+P+y) + ttictP-tQiy+Py) + 6cc$y] . (9) 


Similarly, if the error in ( 7 ) is taken in the form 


C 



C q is found to be given by 

c n = 1 - 5 (d_ 1+ d ia 4 +d 2 8 4 ) 
= 1 - 12 (3+2*8 -to; +8) . 


( 10 ) 


Various alternative modes of utilization of the predictor and corrector 
formulas are available in practice. For example, the predictor can be 
used without employing the corrector at all. On the other hand, if the 
corrector is used, it usually is used i te rat ively, with the predictor 
providing the first guess. Qualitatively, some of the arguments for and 
against the various alternatives are as follows: 
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a) Number of derivative evaluations per step. The "predictor-only" 
mode requires only one evaluation per step . If one correction 
is employed, a second evaluation is usually made after the cor- 
rection to e n hance numerical stability. In general, n cor- 
rections require either n or u»l derivative evaluations, 
depending on whether a final evaluation is or is not carried out. 
Evaluations of complicated derivative functions frequently re- 
quire a predominant portion of the total computer time. 

b) Truncation error. Implementation of the corrector reduces the 
truncation error. (It is a simple exercise to show that 

IcJ < |P n l .) 

c) Numerical stability. With regard to both absolute and relative 
stability, the regions of stability become les6 restrictive as 
the number of correction-evaluation iterations is increased. 
Incidentally, these regions become more restrictive as order is 
increased . 

d) Availability of mesh criteria. More effective procedures for 
automatically selecting the mesh increments can be developed 
for some modes than for others. This consideration favors a 
predictor- corrector mode with at least two applications of the 
corrector . 

An empirical program was carried out whereby the various modes were com- 
pared in the actual numerical solution of selected differential equations . 
The mesh increments were selected in a manner such that the total number 
of derivative evaluations was the same for each mode. This work is not 
reported in detail here since an even more extensive testing program of 
a similar nature for the case of fixed step size was carried out and re- 
ported in detail by Hull and Creemer [ 5 ] • There conclusions, favoring the 
mode p-d-c-d-c, are in agreement with those reached in the present study. 
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(Here p denotes predictor, d derivative evaluation, and c corrector.) 
Consequently, the discussion in the remaining sections will be directed 
primarily toward this mode. 


5. NUMERICAL STABILITY FOR SINGLE DIFFERENTIAL EQUATIONS. The con- 
dition of stability referred to in the previous section does not guarantee 
numerical stability for h > 0. A more appropriate analysis of numerical 
stability is presented here. 


First note that each corrector iteration is performed according to the 
equation 


.(k+1) _ 


n+1 


= y„ 


hd_ 1 F(: 


,(k)> 


n+1’ n+1' 


hS d.y 7 
i=0 i n-i’ 


where the superscript k denotes the iteration. Subtracting this 

equation from (7) and employing the mean value theorem gives 


r n+l 


- «^ 1> - 


-c (k) ) 
n+1 n+1'' 


where 


y=*T 


.00 


for some 77 between y n+1 and c^ + j. Thus the following condition is re- 
quired for convergence of the corrector iterations: 


lAd.jk 1 


( 11 ) 
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It is assumed that condition (ll) is satisfied in the following discussion, 
and in fact this condition will be used in the mesh selection procedures 
described in the next section. 


It is also assumed for the purpose of the nvmerical stability analysis 
that 1 is constant, a standard assumption in the literature for fixed 
step size. By appropriate choice of h at each step, X can be made 
nearly constant in the variable mesh case. In practice, however, this 
assumption is usually violated with fixed mesh methods, even when proce- 
dures to frequently double or halve the step size are Included. Further- 
more when numerical, stability is the controlling factor, it is good policy 
to keep h as large as possible without forcing X beyond its limitation 
imposed by the threat of Instability. Bius in this case, the mesh incre- 
ments vised are actually considerably subopt imal at most steps with fixed 
mesh methods . On the other hand, the variable mesh feature obviously 
allows much better optimization when the integration is stability limited. 
Of course when it is not stability limited, variations in X are 
inconsequential . 


Initially let us consider the mode which employs a prediction and k cor- 
rections with a derivative evaluation after each prediction and correction. 

Let C denote the propagated error, y(x ) - c^). Then it can be shown 
n n n 

that satisfies the difference equation 



2 k 

+ C 

i=l n-ij=l 


X J V7 X d. + i e .X k+1 d k 1 b. 

-1 l i=i n-i -1 i 


> 
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except for the fifth order truncation error. The effect of the predictor 
on the propagated error decreases with increasing k because the factor 
(Ad_j)k multiplies the bj in the above equation. In the limit the cor- 
rector alone determines the error propagation, the equation being given 
by 


W 1 "*-*) - 'n< 1+>i o> - ‘.-1*1 - ' °' < 12 > 

In practice, when the mesh increments are small enough to provide a 

reasonably small truncation error, the corrector iterations beyond the 

second are essentially redundant. Hence the above difference equation 

for the propagated error in the corrector alone is adequately representative, for 

practical purposes, of the error propagation for the recommended mode, p-d-c-d-c* 

If the difference equation (12) has constant coefficients, its solution 
€ n can be expressed in terms of the roots p^ of the polynomial equation 

p 3 (l-Xd_ 1 ) - p 2 (l+Ad o ) - pXd 1 - Xd 2 = 0 (13) 

by € n = + + k^p^ (slightly modified in the case of a multiple 

root), where the are constants. Equation (12) has constant coefficients 

as required provided the d^ are constant as well as X. The d^ are constant 
in the case of fixed mesh. In the variable mesh case, it is this investi- 
gator's experience that the d^ vary very slowly when the integration is 
stability limited. This is due to the fact that the ratio a of mesh in- 
crements from step to step remains nearly constant, and the d^ are con- 
stant when the mesh parameters Of, and y are constant. (When CL is con- 
stant, # and y are the constants cl + a? and CL + 0^ + OC^, respectively*) 

Thus it is reasonable to add the assumption of constant d^, for the pur- 
poses of the stability analysis only, and in view of the above remarks it 
becomes convenient to treat numerical stability in terms of the two parameters 
X and a* 
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When X = 0, the fundamental root of the characteristic equation (13) is 
unity and the others are zero. When X ^ 0, one or both of the latter 
roots may become larger in modulus than the fundamental root. This is 
a condition of relative numerical instability [6], whereas absolute 
numerical instability occurs whenever any root is greater than one in 
modulus or when a root of unit modulus is a multiple root. Applying 
these conditions as definitions, regions of both relative and absolute 
stability have been computed by tracking the roots of (13 ) - These re- 
gions are shown in Fig. 1 in terms of the parameters X and &. Although 
it is interesting to note the behavior for very large and small a, in 
practice actually remains fairly close to unity. Also shown in 
Fig. 1 are the curves X d_^ = ±1, which indicate the region for which 
the corrector iterations converge, and within which the stability regions 
have meaning. 


k. NUMERICAL STABILITY FOR SYSTEMS OF DIFFERENTIAL EQUATIONS • 
The variable mesh formulas are applicable for systems of differential 
equations of the form 


djr 

dx 


(i) 


P ± (x,y (l) ,y <2) , 


y (N) >, i - i. 


2 , 


( 1*0 


In considering numerical stability for this case, equation (12) for the 
propagated error is replaced by 

(i-d^hG)^ - (l+d o hG)c n - - d 2 hG€ n _ 2 = 0, ( 15 ) 


v n' n 

identity matrix and G is the Jacobian matrix with elements G 


where denotes the vector with components y 




I is the 

- SF./dy( J ) 

which are assumed constant, as in the case of a single equation. A cursory 
analysis of numerical stability is available through consideration of a 
characteristic polynomial corresponding to a majorization of equation ( 15 ). 
However, a more detailed approach involving the eigenvalues of the matrix G 
has been pursued in the present study. 
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Pre-multiplying equation (15) by a matrix T, representing a nonsingular 
linear transformation such that TGT = J is in canonical form, gives 

(l-d_lhj)^ n+ i - (l+d o hj)^ n - d 1 hJ^_ 1 - d 2 hJ^_ 2 = 0 , (16) 

where 7 ^ = Te^. The diagonal elements of J are the eigenvalues of G, 
and if these are distinct, all the off diagonal elements of J are zero. 

In this case the system of difference equations for the propagated errors 
becomes uncoupled in passing from (l 5 ) to (l 6 ), and the relevant character- 
istic polynomial equation is again given by (13), with X talcing on the 
values hJ\^ . If the eigenvalues of G are not distinct, the analysis is 
more complicated, as indicated in [ 1 ], but the results are essentially 
the same. In either case, however, Fig. 1 is inadequate because some 
of the J\_. may have nonzero imaginary parts. 

It is easy to show that the zeros of any polynomial whose coefficients 
are themselves polynomials in a complex variable X are the complex con- 
jugates of the zeros of the same polynomial with X replaced by its con- 
jugate. Thus we need only track the roots of ( 13 ) for values of X with 
positive imaginary parts, the regions of numerical stability in the lower 
half of the X^-plane then being given by symmetry. 

The problem of determining regions of stability for fixed a. bas thus 
been reduced to computing the roots of (13) for incremental values of X 
in the upper half ^-plane and deciding at each point whether or not we 
have stability according to some appropriate definition involving the 
roots. \fe will limit ourselves to relative stability. 

Choosing a definition of relative numerical stability presents an in- 
teresting situation. (We ignored this situation in the case of a single 
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differential equation. It was present but rather inconsequential.) One 
would like a definition which not only provides a unique decision regard- 
ing stability at each point but also reflects one’s intuitive notions of 
relative stability. For example, it is distressing to find it possible 
to pass repeatedly back and forth from stability to instability as IX| 
increases along some specified path. Two definitions were considered in 
the present study — one an extension of the Ralston definition used above 
for single differential equations, and the other a definition used by 
Crane and Klopfenstein [7] and also by Krogh [8]. Both definitions lead 
to meaningless relative stability boundaries for fairly large complex X. 
As a practical matter, however, it should be remembered that numerical 
stability is irrelevant for sufficiently large X since either the trunca- 
tion error becomes prohibitively large or convergence of the corrector 
iterations is not obtained. 

The generalization of Ralston’s definition to apply to systems was con- 
sidered by Lea [9]* Lea defined the principal root of the characteristic 
polynomial equation as the continuous function of h satisfying the poly- 
nomial equation and talcing on the value unity at h = 0. All others were 
called extraneous. Actually however, this "definition” fails to distin- 
guish between the principal and extraneous roots because two of them may 
satisfy the requirements of the principal root. The following example 
illustrates this deficiency and further illustrates the inability to de- 
cide between stability and instability for a particular value of X. 

For ql = 1 (fixed step size) the three roots of equation (13) are shown 
in the p- plane (Fig. 2). The values corresponding to X = (-1,2) are in- 
dicated by circles. Moving from the origin in the Xr-plane counterclock- 
wise around the rectangle to (-1,2), the roots proceed in the p- plane 
from the points (l,0), (0,0), and (0,0) to the circled points along the 


15 




Figure 2. Example of Ambiguity In Relative Stability Definition 



paths indicated by the arrows. The point X = (-1,2) appears stable 
according to the Lea definition since the root which started at (l,0) 
is the largest. However, as we continue around the rectangle in Xr-plane 
we see upon returning to the origin that the root which started at (l,0) 
is now at (0,0), while one of the roots which started at (0,0) is now at 
(l,0). In other words, if we had proceeded clockwise in the X-plane, the 
point X = (-1,2) would appear unstable. 

This problem does not develop with small X; that is, when we consider a 
somewhat smaller rectangle the roots return to their starting points. 

On the other hand, the problem does preclude a complete partitioning of 
the X^piane into meaningful regions of stability and instability by this 
procedure . 

The alternate definition does uniquely partition the X~plane into regions 
of stability and instability, but these regions are not acceptable for 
large X . The problem here, although not recognized in either f7] or 
[8], is the one mentioned earlier of alternating between stability and 
instability as X increases. According to this definition, a method is 
relatively stable if the modulus of each of the roots, other than the 
one nearest exp(X), is less than or equal to exp[Re(X)], with equality 
permitted for simple roots only. 

To illustrate the problem with this definition we note first that for 
a = 1 , the roots of equation (13) go from the "aource points,'* (l,0), 

(0,0) and (0,0), to the "sink points," approximately (-2.37,0.0), 
(0.13,-0.17), and (0 .13, +0. 17) , not necessarily respectively, as X goes 
from the origin to infinity along any path in the X^-plane. Consider now, 
for example, X moving along the real axis to (0.5,0.0) and then vertically 
to infinity. For the vertical portion, exp(X) traverses again and again 
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the circle in the p-plane with radius exp(0.5) and center (0,0). Event- 
ually, when the three roots are sufficiently close to their sink points, 
they are each nearest exp(X) for a portion of each cycle of exp(X). Thus 
by definition the method is relatively stable for the portion of each 
cycle when the root near (-2.37,0.0) is the closest to exp(X) and unstable 
otherwise. In this manner, on the vertical line Re(X) = 0.5, ve have 
stability up to Im(X) = 3-0, then instability to about 8.3, stability 
again to about 9.9, etc. 

Since the second definition has the practical advantage that its applica- 
tion is independent of path in the X^-plane, and since the problem just 
noted apparently occurs only for excessively large X ? there is no practical 
difficulty in its usage: one simply ignores stable regions lying "outside" 

unstable regions. 

Consequently the results shown in Fig. 3 were obtained by applying the 
second definition. The two definitions give very similar results for 
small X and reasonable values of O', say l/k < OL < 4 . 

Also shown in Fig. 3 are the curves | | = 1. In a manner analogous 

to the case of a single differential equation, it can be shown that for 
the dominating eigenvalue of the Jacobian matrix of the system, the con- 
dition |Ad_^ | < 1 is necessary for convergence of the corrector iterations. 

5 . CRITERIA FOR SELECTING MESH INCREMENTS. An algorithm for the 
solution of differential equations by variable mesh procedures would be 
incomplete without reasonably sound, general purpose criteria for deciding 
what step size to use at each step of the integration. The main informa- 
tion required for specifying effective criteria was developed in the pre- 
vious sections. In essence, the mesh selection procedure discussed below 
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represents an attempt to choose each step size Just small enough so that 
the following three criteria are satisfied in the numerical solution of 
single differential equations: 

a) The relative truncation error must remain within a prescribed 
tolerance 6. 

b) The condition for convergence of the corrector iterations must 
be satisfied. 

c) The method must possess relative numerical stability. 

Let P n+ ^ and c n+ ^ denote the predicted and final corrected approximations 
for y( x n+ i)* Let H = x n +2 “ x n+l e 8 * e P fl i ze *>° b e used in com- 

puting the solution at x , 0 , and let a, be the new value of ol as determined 

n Tii t 

by the truncation error criterion in a manner described below. (Thus 
from the truncation error criterion we will get H = h/o^..) 


Using the truncation error terms for the predictor and corrector formulas 

5 v 

obtained in Section 2, we can eliminate the factor h_ y n and obtain, through 
fifth order in h, the equation 


5! 


y ( x n+l> - c 


n+1 


/ C n+1 p n+l\ 

\ V®. ) 


where P n and C n are given by equations (9) and (10), respectively. 


We 


want to find a.^ such that the relative error in c n+ g is 6, that is, 


| y <W ~ c n+2 | = 6 | >| * 
In practice we actually set 
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giving 


C n( c n+l~Pn+l) 


TT C , (p -C ) 
n+l v n n' 


1/5 


c im-1 


4 o 


(17) 


If c ■ 0, absolute rather relative truncation error must be 

Qfl 9 

considered, the allowable tolerance depending on the range of the machine. 


Criteria b and c above are combined to produce a single value U c for the 
mesh parameter a at the new step. To this end, we solve H = h/br simul- 
taneously with expressions approximating the boundary of the intersection 
of the regions of relative stability and iteration convergence shown in 
Pig. 1. For this purpose the following expressions have been found to fit 
the boundary data accurately: 


f < 0: 

y 


0<a c <.25: Hf y = -3.2a5 /2 

.25 s a < 1.0: Hf = .17 - 1.096- 

c y c 

1.0 s a <»: Hf = 1.08/a -2, 

c y c 9 

Cl q < .25: not permitted (see Fig. l) 

.25 s a c 5 1.0: Hf y = | [2 + (1 - a c ) 7/4 ] 

1.0 <a c Hfy = 2 + 


An approximation for f^ can be obtained from computations from the com- 
pleted step: 


f ft* 

y 




n+l * n+l ' 


*n+l °n+l 
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^ 1.0, an iterative scheme is used to solve 


For the case f > 0, .25 ^ a 
y c 

for a. : 
c 


a 


(i+D 




8 + 4(1 - a 


(i))7/4 ' 
c ' 


It has been determined (by actual calculations) that with 0^°^ = hf /3, 

f 1 ) ^ 

a} ' is always correct to within two units in the second decimal place. 
Thus ol is computed according to the following simultaneous solutions of 
each of the above equations with the equation H * h/tt^ : 


- 00 < hf £ - .92: ol = (1.08 - hf )/2 

y c y ' 

-.92 < hf ^-.025: a = (.17 +V.03 - 4.36hf )/2.18 

y c y' 

-.025 < hf < 0: a c = (-hf 

0 ^ hf £ .875: a = .25 

y c 

.875 < hf < 8/3: d c = 3hfy/[8 + 4(l - hf^) 7 ^] 

8/3 * hf < »: a = (hf - 2/3)/2. 

y c y 

The new step size H can then be taken as h/a, where a = max(a, ,a ) . 

"t c 

Although this policy has proven satisfactory in practice, it is possible 
that it could produce a new step size which is substantially different 
from the preceding one (but not likely because of the contracting charac- 
ter of the fifth root), and this in turn could result in a subsequent 
loss of accuracy. Therefore the writer recommends the addition of a 
precautionary restriction, such as 2/3 < a < 3/2, using a smaller or 
larger interval depending on the requirements of the particular problem 
being solved. 
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Experience indicates that when even only a moderate degree of accuracy 
is required, the numerical solution of most problems is limited by the 
truncation rather than the stability (or convergence of the iterations) 
criterion. Of course it may be that the truncation criterion is limit- 
ing the step size by detecting numerical Instability of the predictor; 
we know for example that the numerical integration of stiff equations 
is limited by stability. At any rate, when we begin to examine mesh 
criteria for large systems of differential equations, it is especially 
fortuitous that satisfying the truncation error criterion usually pre- 
cludes instability, because in this case the truncation criterion is the 
only one which can be feasibly incorporated into the algorithm. For 
large systems the amount of computing time required to evaluate either 
the Jacobian matrix G or its eigenvalues at each step would usually 
be prohibitive. Of course for certain small systems it may not be pro- 
hibitive, and then the results shown in Fig. 3 can he incorporated in a 
manner analogous to that given above for obtaining a c in the case of a 
single differential equation. This procedure has proved successful for 
selected systems although it did not alter the mesh increments substan- 
tially from those selected by the truncation criterion alone when 
reasonably small values of 5 were used in the latter criterion. 

The mesh selection procedure reconmended for most large systems thus con- 
sists of using only the truncation error criterion. Values of are 

computed from equation ( 17 ) for each component of the system, and then 
a is set equal to the fifth root of the largest of these. 


6 . NUMERICAL TESTING AND COMPARISON WISH OTHER METHODS . The vari- 
able mesh multistep method has been tested by applying it to several 
single differential equations and to several systems of differential 
equations. This testing has given a fairly thorough demonstration of 
the effectiveness and reliability of the algorithm. One system of 
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substantial importance for which the variable mesh approach proved 
especially effective was the problem of heat transfer to a super critical 
fluid with variable physical properties and fully developed turbulent 
flow in a smooth tube [l]. Another system, discussed in [10], was a 
stochastic model of enzymatically controlled cooperative unwinding and 
template replication of biological macromolecules . Due to the compli- 
cated mathematical formulation of these problems, they will not be given 
in detail here. However, several simpler test problems are listed in 
Table 1. 

Most of these test problems were selected because of their inherent poten- 
tial, both in the behavior of the solutions end in the behavior of the 
partial derivatives of the right hand sides with respect to the dependent 
variables, for producing numerical difficulties . Some are particularly 
suited to a variable mesh treatment while others. Nos. 5, 6, 10, and 11, 
can be solved efficiently with constant mesh increments . In the latter 
cases it is important to note that the accuracy obtained by the variable 
mesh method was about the same as that obtained using constant increments 
with the same number of steps . This indicates that the variable mesh pro- 
cedures do not have a degrading effect when they are used unnecessarily. 

Each equation was solved on the IBM System 360 using single precision 

starting values and double precision floating point operations to advance 

the solution. Values of 6, the target relative truncation error, 

-6 -1 

ranging from 10 to 10 were used for each equation. lfce accuracy 
obtained was roughly proportionate to the values of & specified. It 
was noted that the step lengths were limited almost entirely by the trun- 
cation error for the smaller values of 6 with the stability/ convergence 
criterion becoming of increasing importance with Increasing 6. 
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TABLE 1 


PROBLEMS USED TO TEST VARIABLE MESH METHOD 


Problem 

— 

Differential 

Equation(s) 

— 

Integration 

Interval 

Initial 

Value 

Exact 

Solution 

1 

y = -40xy 

-1 £ X £ 1 

exp(-lO) 

exp(l0-20x 2 ) 

2 

y' = (2xy) -1 

1 % £ 10 20 

0 

Vln(x) 

3 

y' = y/x - (l/x)cos(l/x) 

-1 * x * -.01 

sin(l) 

x sin(l/x) 

4 

y / = -exp(x) y 

0 * x * 5 

exp(-l) 

expC-exp(x)] 

5 

/ 

y = -y 

0 * x ^ 10 

1 

exp(-x) 

6 

/ 

y = y 

0 s I MO 

1 

exp(x) 

7 

y' = _y/ z 

OMM 

exp(-l) 

expC-exp(x)] 


z - -z 


1 

exp(-x) 

8 

y' = y(y/z+i) 

IT\ 

VI 

H 

VI 

o 

-exp(-l) 

-exp[x-exp(x)] 


z = y 


exp(-l) 

expt-exp(x)] 

9 

y # = y 2 /z-40z 

-1 £ X * 1 

40 exp(-lO) 

-4 Ox exp(lQ-2Qx 2 ) 


z = y 


exp(-lO) 

exp(l0-20x*) 

10 

y' = -2(y+z) j 

0 M MOO 

0 

-2 exp(-x) ain(x) 


z = y 


1 

exp(-x)[sin(x) + coa(x)] 

11 

y* = -eip(-x)-100z 

0 * x * 1.5 

2 

exp(-x)+ exp(-lOOx) 


z / = -100 z 


1 

exp(-lOOx) 

12 

y. = -z/x 

-1 s x ^ -.01 

cos(l)-sin(l) 

sin(l/x)-(l/x)co8(l/x) 

i 

z = y 

i 


sin(l) 

x 8in(l/x) 




Some of these problems, 1, 9, and 12, were used In comparing the new 
algorithm with other fourth order numerical methods which also permit 
some variability in the mesh increments. The other methods used were 
the standard fourth order Runge-Kutta method, the Nordsieck method, and 
the basic Adams -Bashf orth/Adams -Moulton method, allowing doubling and 
halving of the increments with the latter. As indicated below, the new 
method proved superior to the other methods for these problems. 

Since the Runge-Kutta method requires four derivative evaluations per 
step while the others were used with only two evaluations per step, half 
as many steps were used with the Runge-Kutta method as with the other 
three. For this method the step sizes were obtained by linear interpola- 
tion of an input table, re-applying the method with different tables 
until no improvement could be obtained. 

The Nordsieck method permits Increasing (or decreasing) the step size by 
a factor 6 (or 1/9) . The test problems used in the present study were 
solved with 0 = 2, the value emphasized in [3] where the symbol "0" is 
used for this factor, and also with smaller values to permit more gradual 
varying of the increments. In addition, Nordsieck' s interval control 
mechanism requires a parameter "e" used in a manner to imply a target 
error 0” e . For each value of 0, the problems used here were solved 
with several values of e, seeking one which produced the number of steps 
commensurate with the number used by the other methods. However for 
0 = 2, the Nordsieck method used too many steps even when e was reduced 
to unity. (In fact, considerable difficulty was encountered in trying to 
locate values of 6 which were usable in this sense. Successful choices 
are indicated in Table 2 .) It is also noted here that it was not neces- 
sary to use Nordsieck' s starting procedure for the test problems since 
all the required initial Information was available. 
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For Problem 1, the absolute value of the relative error In the solution 
obtained by each of the four methods is shown in Fig. 4. For this problem, 
the entries in Table 2 are the areas under the curves of Fig. 4. For the 
other two problems, the entries in Table 2 reflect alternative measures 
of relative error which are more appropriate for the numerical solutions 
obtained for those two systems of equations. As can be seen from the 
table, the new variable mesh method gave the best performance; and the 
basic Adams method, augmented with interpolation procedures to permit 
doubling and halving, also did considerably better than the other two 
methods . 


TABLE 2 

COMPARISON OF RELATIVE ERROR 


Problem 

Variable Mesh 

Adams 

Runge-Kutta 

Nordsieck 

1 

1.5 x 10~ 4 

2.3 x 10 -4 

2.8 X 10 -3 

3.0 X 10 -3 (0 = 1.01) 

9 

1.8 x 10 -3 

6.0 x 10 -3 

5.2 x 10 -2 

1.7 x 10 -2 (0 = 1.01) 

12 

2.0 x 10 -2 

5.5 x 10 2 

1.0 x 10 -1 

6.5 x 10 -1 (0 = 1.5) 
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7. EXTENSION OF THE METHOD FOR STIFF EQUATIONS • In the computational 
work carried out with regard to the basic variable mesh method, it was 
observed that the mesh increments are usually limited by the truncation 
error criterion. But numerical stability may be the limiting mechanism in 
the case of large Lipschitz constants. This may arise, for example, in 
studies of mechanical, electrical, chemical, or biological processes where 
a system of differential equations is used to represent the simultaneous 
relaxation of the various components at greatly different rates . 

Several investigators have considered the problem of stability limited inte- 
gration. Crane and Klopfenstin [ 7 ] and Krogh [8] have attacked the problem 
through optimum choice of the coefficients In the "predictor" of standard 
predictor-corrector type methods. They claim slight Improvement for sta- 
bility limited problems, but this is at the expense of substantial degrada- 
tion for other problems . Transformation methods have also been proposed, 
but these are time consuming and difficult to apply in general. A more 
promising approach appears to be that followed by Treanor [ll] where the 
method is designed from the outset with this type of problem in mind. In 
the following paragraphs an algorithm is presented which is similar to that 
of Treanor but has the advantage of relying on the variable mesh, multistep 
framework rather than the Runge-Kutta framework used by Treanor. 

The central idea is to extend the method by including explicit dependence on 
the dependent variable as well as on the independent variable in the local 
approximating function for the derivative, in particular, let us assume 
that f can be approximated in a neighborhood of (x^, y^) as linear in 
y and cubic In x: 

f(x, y) k n » - P(y - y n ) + A + B(x - x q ) + | (x - x n ) 2 + | (x - x Q ) 5 , (18) 
where the coefficients P, A, B, C, D are to be determined. 


29 


£2 



The equation dy/dx = can be integrated in closed form from x q to 

Py 

x &+1 , using the integrating factor e and integration by parts, to 
give a formula for computing y n+1 * The result is 

y n+l " y n + A<1 + (h " ^ (B / P " G/p2 + D / p5) 

+ h 2 (C - D/P)/2P + h 5 D/6P , 
where q = (1 - e" Pll )/p. 

It is the presence of the coefficient P that provides the great advantage 
with regard to numerical stability. If the usual type of analysis of 
numerical stability is applied to this method, numerical stability is in- 
dicated for arbitrarily large |df,/dy|. 

If P is put equal to zero In (18), then the above formula for y n+ ^ is 
replaced by 

y n+l = y n + Ah + Bh 2 /2! + Ch 5 /3! + DhV^! , 

which can be shown to be identical with the basic variable mesh method, 

provided the coefficients are determined appropriately, as follows: The 

five coefficients in (l) axe determined by equating k Q to five previously 

computed values of f . Coefficients for a predictor formula are found by 

matching at the points ( x n-i ^ c & ), i 55 2, 3, and ( x n > P n )j 

where p denotes the predicted value of y(x ) and c denotes the 
n n n 

corrected value . For a corrector, the point ( x n+ p ^n+1^ instead 

° f %-■}'> ■ 
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If the function f is actually at most a linear function of y plus a 
cubic function of x, the formulas give the exact value of y( x n+1 ) • How- 
ever no expression for the truncation error is available in general. 

A computer program was written to test this method with selected differential 
equations. The results of this testing are typified by the following example 
Hie system of equations y' = - 2xy/z, z ' = - 2xz, with initial conditions 
y(0) m z(0) m 1, was solved by both the basic method and the extended 

method in the interval 0 ^ x £ 2 . 2 . This system is not stiff near x «■ 0 , 
but becomes stiff for larger x. The exact solution is y « exp[- exp(x 2 ) ], 

o 

z - exp(- x ). Using the truncation criterion described earlier (with 
& ■ .0001) for mesh selection, accurate solutions were obtained by both 
methods for the dominant component z while the solutions for y remained 
Btable and reasonable (< .0736 error for z and < 10% error for y) . Hie 
basic method required 59k steps compared to 270 for the extended method. 

As expected, all the advantage provided by the extended method came in the 
latter stages of the integration. 

8. APPLICATION TO PARTIAL DIFFERENTIAL EQUATIONS. Some consideration 
has been given to applying the variable mesh multistep approach to partial 
differential equations. To Illustrate one straight-forward way of accom- 
plishing this, consider a method which, when applied to the simple equation 
Uj. “ A u^, reduces to the formula 

Vi,r\j t4 4t[a -i b 2 Vi,j ♦ 4 o D \j 

(19) 

* “l D \-l, J + 4 2 > 
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where the second difference operator D is defined by 

p p 

D u , = (u , - 2u , + u . , )/Ax . 

n, J n,j-l y/ 

With the d^ defined by ( 8 ), (19) can be regarded as an Interpretation of 
( 7 ) • Equation (13) determines the numerical stability, and as can be seen 
from Figure 1 , formula ( 7 ) is absolutely stable for 

- 3 s X s 0, (20) 

for the case a=l, P =* 2, 7 = 3 . 

Now let € . = u(t , x.) - u . . It follows that 

n, J n j n, j 

U ** ( V V - D \ j ' ( 'n,W - 2 *n , 3 + *n,J-l )/4x2 * 

Since (19) is fourth order accurate in time and second In space, we obtain 

2 


* - A At 

n+1 ^ Ax 2 


X d k^ *n-k. 


k = -1 


j+1 ” 2< n-k,j + *n-k, j-1^ * 


Putting c a p n e 1 m we have 

n, j+m 


,3 = p 2 + ^t (e i Ax . 2 + e -i Ax } J 

Ax 2 , 


2-k 


k a -1 


Now if we let 


X = ^4 (e 1 - 2 ♦ o' 1 “) 

Ax 


4 -A At . 2 Ax 


Ax 


2T si n T , 


we deduce from (20) that formula (19) is stable for 
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Ax 4 

for the case a « 1, |B » 2, V-5» 

The above result indicates that this multistep approach has better stability 
than the veil known forward difference method (which requires A At/Ax < l/2), 
allowing a 50 $ Increase in the step size At for given Ax. Furthermore, 
the forward difference method is only first order accurate in time and second 
in space. 

In practice, numerical stability is usually of more concern than local 
accuracy when solving equations of the type considered. Since methods are 
presently available which permit arbitrary mesh ratios, it is felt that a 
more promising alternative to the above is to consider combining the varia- 
ble mesh approach with such unconditionally stable methods . Methods which 
appear to be particularly good candidates for this purpose are the three- 
time- level methods (see e.g. [12 ]) which also permit accurate hn.ndi ing of non- 
linearities in the differential equations. A program along these lines is 
planned. 
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